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[i] Estimation of ice sheet mass balance from satellite altimetry requires interpolation of 
point-scale elevation change (, dH/dt ) data over the area of interest. The largest dH/dt values 
occur over narrow, fast- flowing outlet glaciers, where data coverage of current satellite 
altimetry is poorest. In those areas, straightforward interpolation of data is unlikely to 
reflect the true patterns of dH/dt. Here, four interpolation methods are compared and 
evaluated over Jakobshavn Isbrae, an outlet glacier for which widespread airborne 
validation data are available from NASA’s Airborne Topographic Mapper (ATM). The 
four methods are ordinary kriging (OK), kriging with external drift (KED), where the 
spatial pattern of surface velocity is used as a proxy for that of dH/dt, and their 
spatiotemporal equivalents (ST-OK and ST-KED). KED assumes a linear relationship 
between spatial gradients of velocity and dH/dt, which is confirmed for both negative 
(Pearson’s correlation p < —0.85) and, to a lesser degree, positive (p = 0.73) dH/dt values. 

When compared to ATM data, KED and ST-KED yield more realistic spatial patterns and 
higher thinning rates (over 20 m yr 1 as opposed to 7 m yr 1 for OK). Spatiotemporal 
kriging smooths inter-annual variability and improves interpolation in periods with sparse 
data coverage and we conclude, therefore, that ST-KED produces the best results. Using 
this method increases volume loss estimates from Jakobshavn Isbrae by up to 20% 
compared to those obtained by OK. The proposed interpolation method will improve ice 
sheet mass balance reconstructions from existing and past satellite altimeter data sets, with 
generally poor sampling of outlet glaciers. 
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1. Introduction 

[ 2 ] Three main approaches exist to estimate the mass 
balance of the Greenland and Antarctic ice sheets, each with 
their own errors and uncertainties. First, mass balance may 
be estimated from the difference between the net surface 
mass balance and ice flux through a gate [Rignot and 
Kanagaratnam, 2006; Rignot et al., 2008b]. Secondly, var- 
iations in the gravity field can be used to estimate mass 
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changes of the ice sheets using the GRACE satellites 
(Gravity Recovery and Climate Experiment) [Wouters et al, 
2008; Luthcke et al., 2006]. Thirdly, mass changes can be 
estimated from elevation change (dH/dt) measurements 
using radar or laser altimetry [e.g., Zwally et al., 2005, 2011; 
Sorensen et al, 2011], 

[ 3 ] To obtain an estimate of volume change, point-scale 
dH/dt measurements need to be integrated over the basin and 
thus require spatial interpolation. A method that is usually 
employed to achieve this is ordinary kriging (OK) [e.g., 
Zwally et al., 2005; Bamber et al., 2001], Typically, the 
largest dH/dt values occur in the relatively narrow parts of 
outlet glaciers with high velocities [ Krabill et al., 2004; 
Joughin et al., 2008a; Pritchard et al., 2009]. With con- 
ventional satellite radar altimetry, however, it is difficult to 
obtain sufficient measurements in these areas: because the 
radar beam width is about 20 km [ Thomas et al., 2008], 
narrow outlet glaciers (such as Jakobshavn Isbrae) cannot be 
adequately sampled, and steep slopes cause a displacement 
of the signal, as the return signal does not originate from 
nadir but from the highest point within the radar footprint 
[e.g., Brenner et al., 2007]. While this slope-induced error 
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can be corrected for to some degree, small-scale relief can 
result in noisy height estimates and loss of lock \Bamber and 
Gomez-Dans, 2005], Another issue is the spatial coverage; 
the horizontal resolution is bounded by either crossover 
points (exact locations where ascending and descending 
satellite tracks cross each other) or across-track spacing 
[Thomas et al., 2008], These problems degrade performance 
in the areas where the largest changes occur: the relatively 
narrow outlet glaciers along the ice sheet margin. For future 
studies, the recent launch of CryoSat-2 may partially alle- 
viate these problems due to its higher resolution and latitu- 
dinal coverage (up to 88°N). Satellite laser altimetry can 
give an increased accuracy over outlet glaciers, but the laser 
instrument on the ICESat satellite (GLAS; the Geoscience 
Laser Altimeter System) suffered from operational problems 
[Abshire et al., 2005], resulting in reduced spatial and 
temporal resolution. Stereo photogrammetry can be used 
to determine volume change locally [e.g., Steams and 
Hamilton, 2007], but its accuracy and applicability limits 
it to marginal areas with adequate contrast and large eleva- 
tion differences. It is not suitable for ice sheet wide volume 
change estimation. 

[ 4 ] Thus, given the limited number of valid dH/dt mea- 
surements in the high-velocity parts of the ice sheet (greater 
than 100 m yr -1 ), it is unlikely that an interpolation obtained 
by, for example, OK will reflect the true spatial pattern of 
dH/dt. The similarity of the spatial patterns of surface ice 
velocity and dH/dt [e.g., Rignot et al., 2008a; Joughin et al, 
2008b; Pritchard et al., 2009] suggests that velocity, which 
is available from radar interferometry for nearly all of the 
major outlet glaciers [Joughin et al., 2008a; Howat et al., 
2007], contains information that can aid the interpolation. 
In this paper, we employ an alternative interpolation method, 
kriging with external drift (KED), that takes into account 
auxiliary data [Deutsch and Journel, 1992], KED has been 
applied in geo-sciences in several studies: for example by 
Snepvangers et al. [2003] in interpolation of soil moisture 
values, where net precipitation was used as the external drift, 
or by Schuurmans et al. [2007] for combining precipitation 
measurements from radar with rain gauges. An additional 
extension to the interpolation method is the use of spatio- 
temporal kriging, which handles irregular sampling by 
interpolating both in space and time [Snepvangers et al, 
2003; Kyriakidis and Journel, 1999; De Cesare et al., 
2001]. The resulting interpolation technique, spatiotempo- 
ral kriging with external drift (ST-KED), has two advantages 
compared to OK: first, it makes use of more measurements 
because it samples the past and the future with respect to the 
period of interest; and secondly, velocity gradients constrain 
the interpolation on fast-flowing regions where there are no, 
or limited, data. 

[5] As a case-study area to test whether ST-KED indeed 
improves the estimation of volume change, we use Jakob- 
shavn Isbrse, Greenland’s most active outlet glacier. Jakob- 
shavn Isbrse has been the subject of many investigations 
regarding its recent acceleration and thinning [Joughin et al, 
2004; Thomas et al, 2003; Motyka et al., 2010, 2011; 
Luckman and Murray, 2005], and has been densely sur- 
veyed by airborne laser altimetry [Abdalati and Krabill, 
1999; Krabill et al, 2000, 2004]. The measurements by 
the ATM, mostly located in the high-velocity region, are 
ideal for independently testing our approach. Using these 


validation data, we compare the interpolated patterns of 
dH/dt, thinning rates and errors. 

2. Study Area and Data 

2.1. Study Area: Jakobshavn Isbrae 

[6] Jakobshavn Isbrae is Greenland’s largest outlet glacier, 
draining about 5.4% of the ice sheet (92,080 km 2 ) [Motyka 
et al, 2011], Its location and size are shown in the inset of 
Figure 1. The glacier converges into a 4 km wide, deep (over 
1000 m below sea level) trunk, flowing into a deep fjord 
where, until recently, it fonned a floating ice tongue of about 
15 km length. Between 1850 and 1950, the position of the 
calving front of this tongue retreated some 25 km along the 
fjord. After 1950, the calving front was more or less stable 
until about 1998 [Thomas et al, 2009; Csatho et al., 2008], 
After 1998, rapid melting of the floating ice tongue caused a 
retreat of the calving front, until it finally disintegrated in 
2003 [Joughin et al., 2004]. This collapse was attributed to 
enhanced inflow of warm water into the fjord and subse- 
quent increased basal melting of the tongue [Holland et al., 
2008; Motyka et al, 201 1]. 

[ 7 ] Along with the retreating calving front, measurements 
of surface velocity [Joughin et al., 2004, 2008a] and dH/dt 
[Krabill et al., 2004] indicated acceleration and thinning of 
the high-velocity part of the glacier. While this area had been 
thickening prior to 1997 [Abdalati et al., 2001; Thomas 
et al., 2003], it has been thinning since then, with rates in 
excess of 1 5 m yr - 1 near the grounding line. Over the same 
period, Joughin et al. [2004] measured ice velocity, aver- 
aged over some 20 locations in the high-speed area of 
Jakobshavn Isbrae. Between 1985 and 1992, velocities were 
slowing down from 6.7 km yr -1 to 5.7 km yr -1 . After 1999, 
however, the glacier has been speeding up, to 9.4 km yr -1 in 
2000 and 12.6 km yr -1 in 2003 [Joughin et al., 2004]. 

2.2. Data 

[s] As input for our interpolation, we employ data from 
two spacebome altimeters: i) the RA-1 (Radar Altimeter) on 
ERS-2 (European Remote Sensing satellite); and ii) the laser 
altimeter (GLAS) on ICESat. In addition, we use data from 
an airborne laser altimeter (NASA’s ATM) for validation 
purposes. ESA’s ERS-2 was launched in 1995 and has been 
in an orbit at an altitude of 780 kilometer with a 35-day 
repeat cycle until September 2011 when it was decommis- 
sioned. In this study, ERS-2 data between 1995 and 2003 
have been used. dH/dt values derived from ERS-2 are based 
on crossover points that were averaged into clusters [e.g., 
Li and Davis, 2008], shown as blue triangles in Figure 1. 
ICESat was launched in 2003 and carried a laser altimeter 
system (GLAS). Due to rapid laser degradation, the system 
was switched on for only two to three campaigns of 33 days 
per year [Abshire et al., 2005], The spacing between tracks 
is therefore relatively large, about 20 km for Jakobshavn 
Isbrae (Figure 1). Along-track spacing is approximately 
1 72 m, with a footprint size of about 60 meters [Zwallv et al, 
2002]. We use ICESat data (release 31) from 14 epochs: from 
February/March 2003 until February /March 2008 [Zwallv 
et al., 2010], The same preprocessing was carried out as 
in Sorensen et al. [201 1]. ICESat repeat tracks are shown in 
red in Figure 1. ATM is also a laser altimeter, but airborne. 
Since 1993, repeated flights have been performed over 
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Figure 1. Overview of the study area and spatial coverage of the data used in this study, overlaid on a 
MODIS image of the area from June 2002. Ice velocity is shown as contours with 200 m yr -1 intervals. 
The thick black line shows the approximate terminus position in June 2004 (from Csatho et al. [2008]). 
Only locations for which a valid dH/dt estimate (see text) could be obtained are shown. The inset shows 
location and size of Jakobshavn Isbras (black shade) and the red square is the extent of the figure. 


Greenland [Krabill et al., 2000, 2004], with sampling over 
parts of Jakobshavn Isbras [Abdalati and Krabill, 1999] 
shown in Figure 1. ATM measures ice-surface elevation 
with an accuracy of 10 cm or better [ Krabill et al., 2000], 
within a swath of 150-200 m at a typical aircraft altitude of 
500 m [ Thomas et al., 2009]. Here, a processed version of 
the data was used, in which footprints (with size 1-3 m) 
were resampled to 70 m overlapping “platelets”, one on 
each side of the aircraft [ Thomas et al, 2009]. 

[ 9 ] Ice velocity is derived from a combination of inter- 
ferometric synthetic aperture radar (InSAR) and speckle 
tracking [Joughin, 2002], In the winter of 2000/2001, 
Canada’s RADARSAT-1 acquired near-complete coverage 
of the Greenland ice sheet. From 2005/2006 this coverage 
was repeated almost annually [ Joughin et al., 2010], We use 
a mosaic of RADARSAT-1 derived velocities averaged over 
the years 2000, and 2005-2008. To define the outline of the 
drainage basin of Jakobshavn Is bra;, finally, the delineation 
of Rignot and Kanagaratnam [2006] is followed. 

3. Methodology 

3.1. Obtaining dH/dt Values 

[ 10 ] To be able to combine all data sets described before, 
we define a grid in polar-stereographic projection (central 
latitude 71°N; central longitude 39°W), at a 1 km resolution. 
The processing approach for laser altimetry (ATM, ICESat) 
differs from that for radar altimetry. The radar altimetry data 
were already preprocessed as crossover clusters, but not 


corrected for the slope-induced error, i.e., the fact that the 
return signal does not come from the point at nadir, but from 
the point within the rather large footprint that is closest to the 
antenna [e.g., Brenner et al., 1983], Near the ice sheet 
margin, the difference between these points can be consid- 
erable (about 14 km for a 1 degree slope). This is corrected 
for using slope and aspect from a digital elevation model 
(DEM) [Bomber et al., 2001] as described in Hurkmans 
et al. [2012], 

[ 11 ] For ICESat, several approaches have been used to 
extract dH/dt [Howat et al., 2008; Pritchard et al., 2009; 
Slobbe et al., 2008]. Because ICESat tracks almost never 
exactly repeat, a regression approach can be used in which 
slope (both across-track and along-track) and dH/dt are 
simultaneously solved for. To estimate dH/dt from ATM, 
footprints are resampled into 70 m platelets as mentioned in 
Section 2.2, and elevations at platelet centers of overlapping 
flightlines are compared [Krabill et al., 2004; Thomas et al, 
2009]. However, in all these studies only one average dH/dt 
was obtained, where temporal resolution was sacrificed in a 
way to increase spatial coverage [Pritchard et al., 2009]. 
Because we are interested in the change of thinning or 
thickening rates through time, we need a higher temporal 
resolution as well. 

[ 12 ] First, for each repeat track (or flightline for ATM), we 
average all footprint elevations within a given 1 km grid cell. 
We again use slope and aspect from a DEM [B amber et al., 
2001] to project all elevation footprints (for ICESat) or 
platelet centers (for ATM) toward the center of the pixel. 
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Figure 2. The dH/dt from all sources plotted over the same MODIS image as in Figure 1 . Elevation rates 
are determined using (a) all data after 1999, or (b-f) data from individual 3-year periods centered on 1998, 
2002, 2003, 2004, or 2006 respectively. Also velocity is shown in black 200 m yr^ 1 contours (up to 
1,400 m yr -1 ). For clarity, only the area around the high-velocity part of the glacier is shown, similar to 
the area shown in Figure 1. The blue star denotes the measurement with the highest thinning rate in the 
image. 


This approach is likely to be less accurate than a method that 
estimates slope directly from the data as was done in, for 
example, Pritchard et al. [2009]. However, we expect this to 
add noise to the results rather than a bias. Because we 
employ rather strict data culling criteria (see below), we do 
not believe our results are significantly influenced by this. 
Time-averaged dH/dt values are also calculated for each 
satellite (using regression including all available data). For 
ICESat, the spatial pattern of overall dH/dt for the entire 
Greenland ice sheet (2003-2007) compares very well with 
those reported by Pritchard et al. [2009] and Sorensen et al. 
[2011] (not shown). 

[ 13 ] For both radar and laser altimetry, we then calculate 
annual values for dH/dt at each grid cell for every year 


between 1997 and 2008 using linear regression. For each 
year, data are used from a three-year period centered on the 
year under consideration, because the temporal sampling, 
especially before 2003, is too poor to allow reliable dH/dt 
calculations within single years. Regression is only carried 
out in grid cells where i) data from at least five overpasses 
are available that ii) span at least one year and iii) have a 
standard error on the resulting dH/dt rate of less than 
0.40 m yr - 1 . The standard error on the regression coefficient 
SE CO ef (in this case dH/dt rate) is calculated by: 


SP CO ef 


' £*?/(” - 2 ) 


(i) 
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Figure 3. As Figure 2, but now standard errors of the dH/dt displayed in Figure 2 are shown. 


where e is the vector of residuals, n is the sample size, and 
x is the input with mean x. It should be noted that this 
standard error is not equivalent to measurement error, but 
takes into account sample size, and variance of both input 
and residuals of the regression as well. The threshold is 
selected by trial and error to avoid a noisy spatial pattern of 
points that are close together and opposite in sign, usually 
because the regression was based on a small subset of 
overpasses. In addition, prior to regression, measurements 
that are outside the range of average plus or minus two 
standard deviations (based on the entire available time 
series) are removed. The number of data points that was 
removed in the culling process is variable between data sets 
and time periods. For ICESat, typically about 25% of data 
points are removed, mostly because less than five overpasses 
were available or the spanned time range was too short. 
Figure 2 shows estimated dH/dt for different periods, based 
on linear regression. Figure 3 shows the associated standard 
errors of the dH/dt. 


3.2. Kriging Theory 

[ 14 ] A widely used interpolation method in geosciences is 
kriging, sometimes referred to as “optimal interpolation”. 
It is a form of linear regression minimizing an estimation 
variance [Deutsch and Journel, 1992], of which several 
varieties exist. In the following, the basic theory and the 
differences between the main types are briefly explained. 
For more information, see for example Deutsch and Journel 
[1992], Herzfeld et al. [2000], and Hengl et al. [2003], 

[ 15 ] In kriging, a value Z(x 0 ) at an unsampled location can 
be estimated from N samples located at x,- using: 

N 

Z(x 0 ) - m(x 0 ) = ^ - «(*/)] ( 2 ) 

i=l 

where Z is treated as a random field with a trend component, 
m, and a residual component R — Z — m, and \ is the weight 
assigned to each sample location. R is treated as a random 
function with a mean of zero and a covariance C(h). C(h ) can 
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be estimated from the data in the form of a semi-variogram 
1(h): 


= i E i^o-zfe )! 2 o) 

where N(h) is the number of data pairs at distance h, and 
Z(x,) and Z{xj) are the data values at points x- t and Xj, 
which are separated by a distance h. The semi-variogram thus 
indicates the expected squared increment of the data value 
over distance h, and is characterized by three parameters 
[e.g., Herzfeld et al., 2000]: the nugget (point-scale vari- 
ance), the range (largest distance at which there is spatial 
dependence), and the sill (the semi-variance at distances 
beyond the range). The covariance C(h) can then be obtained 
by C(h) = C(0) — 7 (/;) where C(0) is the sill. The sill and 
range parameters are derived by fitting a model to the semi- 
variogram calculated from the data. Typically, an exponen- 
tial, gaussian, or spherical equation is used for this [Deutsch 
and Journel, 1992]. 

[16] In most geostatistical applications, data to be inter- 
polated do not only vary in space, but also in time. There- 
fore, several models have been developed for space-time 
variability [Kyriakidis and Journel, 1999]. One of these is 
the product-sum method [De Cesar e et al., 2001], which is 
more flexible compared to other models, does not require an 
arbitrary space-time metric, and it is relatively straightfor- 
ward to fit models to the sample semi-variogram [Gething, 
2006]. In case of spatiotemporal kriging, a two-dimensional 
sample semi-variogram is obtained using a distance lag h s 
and a time lag h, as dimensions, so that for every h t a spatial 
semi-variogram is obtained and vice versa. A semi-variogram 
model is then fitted to both the spatial (average of all /?,’ s) 
and temporal (average of all h s ’ s) sample semi-variogram. 
Using the product-sum model, the modeled spatial and 
temporal semi-variograms are then combined: 

l(h s ,h t ) = (*iC,(0) + k 3 )y t (h,) + (hC,(0) + k 2 )y s {h s ) 

— k\l s (h s )lt{h t ) (4) 

where 7 s (hj) and 7 ,(h t ) are the spatial and temporal semi- 
variogram, respectively; C,(0) and C,(0) are the spatial and 
temporal sills, respectively; and k\, k 2 , and k 2 are given by: 

( (C v (0) + C,(0) — C sl (0) 

1 C,(0)C,(0) 

, C„(0) - C,( 0) 

<k2 ~ C,( 0) (5) 

, C*(0) - C,(0) 

l 3 C t ( 0) 

where C st (0) is the spatiotemporal sill, which is normally 
taken as the maximum value of the sample space-time semi- 
variogram surface \De Cesare et al., 2002]. The space-time 
covariance is then used for three-dimensional interpolation 
(x, y, and t), resulting in an interpolated surface at every 
pixel and at every desired time step, where also samples at 
earlier or later time steps are used in the estimate. 

[17] The treatment of the trend m is the main difference 
between the various members of the kriging family. In this 
study, ordinary kriging (OK) and kriging with external drift 


(KED) are employed. The system of equations that needs to 
be solved for interpolation is derived by expressing the 
variance of the estimation error in terms of the covariances 
C(/j), and then minimizing it by equating the partial deri- 
vatives with respect to A,s to zero (see Deutsch and Journel 
[1992] for details). In the simplest case, where m is assumed 
uniform, this yields: 

N 

E %C(xi - Xj) = C(xi - x 0 ) , i = 1 , 2, . . . , N (6) 

7=1 

where C(x, — xj) is the covariance function for lag h = x t — Xj 
and C(Xj — x 0 ) is the covariance at lag ( x,- — x 0 ). This case is 
called simple kriging (SK). When other types of kriging are 
employed the principle remains, but the form of the equation 
system in equation (6) changes slightly. In OK, the most 
widely used variant, m is assumed to be constant only for the 
local neighborhood of each estimation point (m(x,) = m{xj)). 
To be able to estimate m locally, equation (6) is extended 
with a constraint requiring the sum of weights to be one: 

' N 

E *)C R (xi-xj) + fi(x 0 ) = C R (xi-x 0 ), 

7=1 

< for i= 1,2,. ..,N (7) 

N 

E'V = 1 

7=1 

where n(xo) is a Lagrange parameter associated with the 
constraint on the weights. Alternatively, a spatial trend can 
be assumed to exist in m. In kriging with a trend (KT), 
sometimes referred to as “universal kriging”, the trend is a 
linear or higher order function of the (x,y)-coordinates. 
In case of kriging with external drift (KED), a relation is 
assumed between Z and a secondary variable Y, so that the 
trend can be calculated from the secondary variable. Y is 
assumed to reflect the spatial variability of the variable to be 
interpolated and to be known at every location in the inter- 
polation grid. In the case of KED, the trend is assumed to be 
linear (only a first-order trend): 

m(x 0 ) = a 0 + a\Y(xj) (8) 


[is] The system of equations to be solved is further 
extended, allowing local and implicit estimation ofr/o and ay. 


< 


N 

E A/C(x,- - xj) + /i 0 (xo) + Ah (xo)Y(xi) 
7=1 

= C(xi — .To) , for i = 1 , . . . , N 

N 


E A 7=1 


7=1 


N 


E A 7 7 ( X >) = 

7=1 


(9) 


In this study, the FORTRAN programs GAMV and KT3D 
from the GSL1B geostatistical package [Deutsch and 
Journel, 1992] were used for calculation of the sample 
semi-variogram (GAMV) and the actual kriging (both OK 
and KED). An extension to GSL1B to allow for spatiotem- 
poral kriging was provided by De Cesare et al. [2002], which 
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Figure 4. Scatterplots of dH/dt versus velocity, for (clock- 
wise from upper left) ATM, ICESat, ICESat+ERS-2, and all 
sensors. For each sensor, data from the entire basin 
(Figure 1) were used, and dH/dt values are based on the 
period 1999-2008 (2003-2008 for ICESat, 1999-2003 for 
ERS-2). Lines and correlation coefficients are based on lin- 
ear least squares regression. 

we also employed here. It should be noted that all kriging 
types (OK, KED, KT, etc.) can be used in both spatial-only 
and spatiotemporal contexts. For example, Snepvangers etal. 
[2003] used space-time ordinary kriging (ST-OK) and space- 
time kriging with external drift (ST-KED) to model soil 
moisture contents. In the remainder of this study, we employ 
and compare four kriging types: OK, KED, ST-OK, and 
ST-KED. 

3.3. Interpolation of dH/dt Values 

[ 19 ] Many studies involving spatial interpolation of satel- 
lite measurements have used ordinary kriging (OK) [Zwallv 
et ah, 2005; Bamber et al., 2001; Herzfeld et ah, 2000], 
Considering the data coverage in Figures 1 and 2 however, it 
is unlikely that ordinary kriging, based on satellite altimetry 
alone, will yield a representative spatial pattern. To illustrate 
this, the maximum thinning rate obtained from ATM data 
(over the period 1999-2008) was about 17 m yr _1 , at the 
location shown by the star in Figure 2. Note that this point is 
not located in the center of the main trunk, and some kilo- 
meters away from the grounding zone, so even higher values 
are possible. The maximum thinning rate as derived from 
ICESat, on the other hand, was only about 7 m yr~' and 
values obtained by radar altimetry were even lower than that. 
A possible alternative interpolation method would be the 
hypsometric approach, in which a relationship is assumed 
between elevation and dH/dt [Arendt et ah, 2002]. However, 
since dynamically induced dH/dt and elevation are not 
clearly related, the hypsometric approach is not further dis- 
cussed here. 

[ 20 ] Ice velocity measurements are available at nearly all 
locations on the ice sheet, especially over fast-flowing outlet 


glaciers. From visual inspection of Figure 2, the spatial 
patterns of dH/dt and ice velocity seem very similar. Here, 
we investigate whether there is information in the ice 
velocity pattern that can aid the interpolation of dH/dt. As 
explained earlier, kriging with external drift (KED) uses 
auxiliary information to determine the underlying trend 
model for the interpolation. The secondary variable (veloc- 
ity) is assumed to reflect the spatial trend in dH/dt up to a 
linear rescaling of units [Deutsch and Journel, 1992]. Spatial 
gradients in velocity should thus be related to spatial gra- 
dients in dH/dt. 

[ 21 ] A spatial correlation between (dynamically induced) 
dH/dt and velocity has been noted before [Joughin et ah, 
2008b; Rignot et ah, 2008a; Pritchard et ah, 2009], In 
fast-flowing areas, that are flowing faster than required to 
maintain mass balance, ice stretches longitudinally, which 
causes it to thin vertically [Rignot et ah, 2008a]. The thin- 
ning increases surface slopes, which further increases 
velocity [Joughin et ah, 2008b]. On outlet glaciers in 
both Antarctica and Greenland, significant thinning (and 
associated speed-up) is confined to the fast-flowing areas 
[Pritchard et ah, 2009]. This is also the case for Jakobshavn 
Isbras [Joughin et ah, 2008b]. Therefore it is reasonable to 
assume that there is a relation between the spatial gradients 
of velocity and dH/dt. Such a relation would also depend on 
glacier width and thickness. Because Jakobshavn Isbrae is 
flowing through a deep trough and its thickness is >1000 m, 
and changes in width are also relatively small, changes in 
velocity are dominant in this respect. In addition, it should 
be stressed that such a relation is not fixed: its parameters 
(f?! and a 0 in equation (8)) are determined implicitly from the 
kriging system of equations (equation (9)) and will be vari- 
able in space and time. Because of this flexibility, different 
spatial gradients of dH/dt can result from the same spatial 
gradient in velocity, constrained by the altimetry measure- 
ments. This allows us to use an average, static velocity field 
in order to maximize spatial coverage. 

[ 22 ] dH/dt is not only caused by ice dynamics but also by 
changes in surface mass balance, which will cause noise in 
any relationship between velocity and dynamically induced 
dH/dt. In the high-velocity area of Jakobshavn Isbrae, how- 
ever, thinning rates are >15 m yr _I , whereas maximum 
SMB values are about —4 m ice yr _1 [e.g., Thomas et ah, 
2003]. Dynamic changes are therefore much larger than 
SMB-induced changes and we can assume that ice dynamics 
are the dominant signal. 

[ 23 ] To assess the correlation between velocity and dH/dt, 
the two are plotted against each other in Figure 4. The plots 
and the highly negative Pearson’s correlation coefficients 
(<— 0.85) suggest a near-linear relationship, from which we 
conclude that velocity can be a suitable secondary variable 
for KED. Figure 4 shows correlation coefficients for the 
period 1999-2008 (2003-2008 for ICESat), but they were 
also calculated for individual years (Table 1). It was, how- 
ever, not always possible to obtain a robust correlation 
coefficient for each year due to paucity of data. After 1999, 
correlation coefficients are reasonably high for all sensors. 

[ 24 ] Interestingly, positive correlation coefficients show 
up in the early nineties, in the period that Jakobshavn Isbrae 
was thickening [Abdalati et ah, 2001] suggesting a dynamic 
origin. Apparently, the areas with the highest velocity after 
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Table 1 . Correlation Between Velocity and dH/dt From All 
Instruments 3 



ATM 

ICESat 

ERS-2 

All 

P 

N 

P 

N 

P 

N 

P 

N 

1994 

0.73 

64 

_ 

_ 

_ 

_ 

0.73 

64 

1995 

- 

- 

- 

- 

-0.23 

77 

-0.23 

77 

1996 

0.48 

38 

- 

- 

0.06 

113 

0.49 

151 

1997 

—0.02 

67 

- 

- 

0.17 

115 

-0.19 

182 

1998 

-0.25 

112 

- 

- 

-0.41 

116 

-0.28 

228 

1999 

- 

- 

- 

- 

-0.37 

115 

-0.37 

115 

2000 

- 

- 

- 

- 

-0.32 

114 

0.04 

123 

2001 

0.26 

10 

- 

- 

-0.32 

114 

-0.07 

124 

2002 

-0.63 

26 

- 

- 

-0.31 

114 

-0.81 

140 

2003 

- 

- 

-0.47 

30 

-0.66 

105 

-0.57 

135 

2004 

-0.82 

34 

-0.53 

4674 

- 

- 

-0.82 

4708 

2005 

-0.81 

33 

-0.79 

6181 

- 

- 

-0.89 

6214 

2006 

-0.81 

33 

-0.84 

6007 

- 

- 

-0.91 

6040 

2007 

- 

- 

-0.80 

3986 

- 

- 

-0.80 

3986 

ALL 

-0.85 

524 

-0.88 

6727 

-0.56 

117 

-0.91 

7368 


Correlation coefficients between velocity and dH/dt (p) and the number 
of measurements (N), for ATM, ICESat, ERS-2, and all data combined. 
Values are shown for individual years, based on data from moving 3 -year 
periods, and averaged over the period 1999-2008 (ALL). Data before 
1999 are excluded to limit the influence of opposing signs due to the 
transition from thickening to thinning in 1998 [ Csatho et al . , 2008; 
Thomas et al., 2003]. 

2000, are also the areas that thickened most before 1997. 
Figure 5 shows the corresponding scatterplots for ATM 
measurements in 1994 and 1996. For both years, essentially 
only one repeated flight line was available, so the number of 
samples is limited, and the significance of the trend is low. 
The presence of some warm summers in the mid-nineties 
[: Thomas et al., 2003] and a simultaneous slow-down of the 
glacier [Joughin et al., 2004] suggest that the thickening was 
dynamic, and these positive correlation coefficients seem to 
confirm that. Between 1997 and 2003 correlations are lower 
because the glacier was in a transition between thickening 
and thinning, and very few valid dH/dt estimates are present 
before 2003. Because a 3 -year window was required to 
calculate dH/dt, the period at the end of the ERS-2 era and 
the beginning of the ICESat era has a limited number of 
measurements. After 2003, when many ICESat epochs can 
be used, there is a strong anti-correlation ( p < —0.80). 

[ 25 ] To assess the value of including velocities in the 
interpolation, both OK and KED were used to interpolate the 
dH/dt data obtained in Section 3.1. In addition, to overcome 
the irregular temporal sampling by the different sensors, 
ST-OK and ST-KED were also used. One approach to do 
this would be to interpolate elevations prior to calculating 
dH/dt. Due to the very different spatial coverage of ERS-2 
compared to ICESat, many areas have coverage for only part 
of the period, causing artifacts in dH/dt around these transi- 
tion periods and noisy results. We therefore chose to calcu- 
late dH/dt prior to interpolation as described in Section 3.1 
and interpolate the resulting annual values. Besides, as there 
is no clear relation between elevation itself and ice velocity, 
incorporating the external drift is not straightforward or 
sensible if elevations are interpolated directly. 

[ 26 ] As a first step of the interpolation process, semi- 
variogram models were estimated (see Section 3.2) using the 
GAMV routine in the GSLIB package (for OK/KED) or the 
extension to GSLIB by De Cesare et al. [2002], after which 


range and sill parameters were estimated using standard 
curve fitting routines. Figure 6 shows all resulting sample 
semi-variograms and models. Corresponding parameters and 
variogram model types (gaussian, spherical or exponential; 
see for example Herzfeld et al. [2000] for details) are shown 
in Table 2. 

[ 27 ] For all kriging methods, we use the sample semi- 
variogram calculated from the satellite data only (ICESat 
and ERS-2), as we intend to use ATM data for validation 
only. The resulting semi-variogram and model are shown in 
Figure 6a. For comparison, the semi-variogram resulting 
from ATM data only is shown in Figure 6b. Because most of 
the ATM data are concentrated in the high-velocity area of 
the glacier (roughly a 50 x 50 km area), semi-variance for 
ATM is only considered and shown for spatial lags up to 
about 50 km. Because ice flow is mainly east-west, and most 
data are organized in (almost) north-south oriented tracks, 
dH/dt is expected to be anisotropic. Therefore, specific var- 
iograms were calculated for east-west and north-south 
directions, shown as gray and orange lines in Figures 6a 
and 6b. However, for lags up to about 60 km (all of the 
ATM data), semi-variograms for specific directions are very 
similar. Only for larger spatial lags do the east-west and 
overall semi-variograms diverge, but not significantly. 
Because distances in east-west direction within the basin are 
much larger than north-south, the north-south variogram 
should only, like ATM, be considered at smaller lags. Based 
on these results, we use the overall semi-variograms in the 
remainder of this study. Moreover, when the anisotropic 
semi-variograms were used for kriging the resulting anisot- 
ropy in the interpolated pattern was greatly exaggerated (not 
shown). 

[ 28 ] Figures 6c and 6d show the spatiotemporal semi- 
variograms. From the sample semi-variogram surface (black 


ATM 1 994 ATM 1 996 



Figure 5. The dH/dt versus velocity, specifically for ATM 
data in (a) 1994 and (b) 1996, when Jakobshavn Isbras was 
thickening. Lines and correlation coefficients are based on 
linear least squares regression. Although few ATM samples 
are available and the signal-to-noise ratio is relatively low, 
a clear relation is present between (dynamic) thickening 
and ice velocity. 
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Figure 6. Semi-variograms and their fitted models, based on (a, c, d) input satellite data and (b) ATM 
validation data. Figure 6a shows the variogram for spatial-only kriging; Figure 6c shows the marginal spa- 
tial (lower set of curves) and temporal (upper set) semi-variograms; Figure 6d shows the semi-variogram 
surface from data (black) and model(red). In all plots, black/gray (black for isotropic, gray for anisotropic) 
lines reflect data, whereas red/orange lines indicate fitted models. All parameter values are provided in 
Table 2. 


surface in Figure 6d), the marginal space and time semi- 
variograms are derived, and models are fitted to them 
(Figure 6c). These are combined using the product-sum model 
[De Cesare et al., 2001] and the red square, representing the 


modeled semi-variogram surface (equation (4)) is obtained. 
The marginal temporal variogram is not particularly struc- 
tured and not a very good fit is obtained. It should be noted, 
though, that the temporal resolution was probably not 


Table 2. Semi-variogram Parameters for the Semi-variograms Shown in Figure 6 a 



Satellite Input 


ATM 


Cst 

Space Marginal 

Time Marginal 

Type 

Range 

(km) 

Sill 

(m 2 ) 

Type 

Range 

(km) 

Sill 

(m 2 ) 

Type 

Range 

(km) 

Sill 

(m 2 ) 

Type 

Range 

(km) 

Sill 

(m 2 ) 

1-D 

Exp. 

497.5 

0.216 

Gau. 

63.5 

7.712 

0.16 

Exp. 

54.4 

0.076 

Sph. 

10.2 

0.071 

E-W 

Exp. 

834.1 

0.448 

Gau. 

86.8 

8.352 

0.22 

Exp. 

210.5 

0.123 

Sph. 

9.0 

0.088 

N-S 

Sph. 

117.4 

0.039 

Gau. 

51.9 

7.809 

0.13 

Sph. 

135.7 

0.048 

Sph. 

12.0 

0.041 


“Parameters (Range and Sill) as well as type (“Exp.”, “Gau.” and “Sph.”, are, respectively, exponential, gaussian and spherical), are shown for the 
isotropic case (1-D) as well as two directions (east-west, E-W and north-south, N-S) in the anisotropic case. 
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Distance along centerline [km] 

Figure 7. (a) Location of the centerline of the main trunk 
(approximate line of highest velocity) of Jakobshavn Isbras 
as well as locations of input (brown dots) and ATM (green 
stars) measurements, (b) The dH/dt versus distance along 
this centerline for all interpolation methods and data points, 
as well as velocity (black line). 


sufficient to capture a complete covariance structure 
because of the regression in moving three-year periods. 

4. Results 

[ 29 ] Based on the parameters from Table 2, annual dH/dt 
estimates for the entire Jakobshavn Isbra area were derived, 
using all four kriging variants. In addition, longer-term 
dH/dt values were interpolated using all data after 1999, i.e., 
1999-2003 for ERS-2 and 2003-2008 for ICESat. Data 
before 1999 were excluded here because the change from 
thickening to thinning took place in around 1998 and 
assuming a linear change rate would thus be invalid. 

[ 30 ] Figure 7 shows temporally averaged dH/dt along a 
transect following the approximate line of maximum veloc- 
ity from the grounding zone toward the accumulation zone. 
dH/dt values are shown for all four kriging methods, as are 
input (ERS-2 and ICESat) and validation (ATM) measure- 
ments along the transect. It should be noted that OK and 
KED use the temporally averaged dH/dt values (which are 
also shown) as input, whereas for ST-OK/ST-KED interpo- 
lation results were averaged. In Figure 7b, therefore, the 
transects for ST-OK/ST-KED do not pass exactly through 
the input data. In addition. Figure 8 shows the spatial 


patterns of dH/dt for the same period. Input data are shown 
as well (Figure 8a), and Figure 8b shows interpolated ATM 
data (using ordinary kriging), which can be considered as the 
“true” spatial pattern because of the high sampling density in 
the area of interest. To get an estimate from ST-OK and ST- 
KED for the same period, the average dH/dt over the years 
1999-2008 was calculated (Figures 8e and 8f). 

[ 31 ] Both Figures 7 and 8 show that interpolated thinning 
rates are never higher than measured ones (about 7 m yr -1 ) 
if no velocity data are used. This is clearly unrealistic 
(Figure 8a). As was mentioned in Section 3.3, ATM data 
show much higher thinning rates. The highest measure- 
ment, close to the grounding line in 2004 (Figure 2e), is 
about 25 m yr -1 , which is an indication of the maximum 
thinning rate near the grounding zone. 

[ 32 ] The interpolated spatial patterns using KED and 
ST-KED are evidently more realistic than those using OK/ 
ST-OK compared to the interpolated ATM data (Figure 8). 
Also the interpolated thinning rates along the transect 
closely follow those of the ATM validation data (Figure 7). 
It should be noted that for the ATM data shown in Figure 7 
the threshold on the standard error was slightly relaxed to be 
able to show more ATM data points (0.6 m yr -1 whereas it 
is 0.4 m yr - 1 in the rest of this study). 

[ 33 ] The velocity along the transect, which is also 
shown in Figure 7, rapidly increases from about 4 km yr -1 
to 12 km yr -1 in the 10 km nearest to the grounding zone. 
ATM measurements along the transect show that the inter- 
polated dH/dt values for KED and ST-KED are accurate 
upstream of that point. Closer to the grounding zone, how- 
ever, thinning rates as predicted by KED/ST-KED follow a 
similar gradient to velocity (as would be expected from the 
external drift algorithm), rapidly reaching improbable values 
(up to —80 m yr -1 ). Two explanations are possible for this: 
first, while we excluded the area west of the approximate 
grounding line in 2004, the grounding line has been 
retreating in recent years [Csatho et al, 2008; Thomas et al, 
2003] and might have migrated several kilometers further 
inland by 2008. Part of the area with very high velocities 
may thus have become ungrounded. Second, the sudden 
increase in velocity gradient suggests a change in flow 
regime. In both cases, the relation between velocity and 
dH/dt breaks down and we thus decided to exclude the area 
with velocities >4000 m yr -1 , and the area within 5 km 
east of the 2004 grounding line. The combined area that is 
removed is relatively small (Figure 8). 

[ 34 ] Statistics of dH/dt averaged over the same period are 
shown in Table 3 . As we are focusing here on the high- 
velocity area, statistics were calculated using pixels associ- 
ated with velocities >100 m yr -1 only (Figure 8). However, 
for the entire basin, average dH/dt values are also negative 
(not shown). This is in line with the net mass loss reported 
by other studies [e.g., Rignot et al., 2008b; Sorensen et al., 
2011 ]. 

[ 35 ] To quantify the errors of all four interpolation meth- 
ods, a cross-validation approach was carried out [Deutsch 
and Journel, 1992], Using a separate run for each interpo- 
lation algorithm, every available altimetry measurement is 
left out and estimated from its neighbors. The error is then 
the difference between the interpolated and the actually 
measured value. Statistics of the errors, for individual 
years and the period 1999-2008, are shown in the form of 
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Figure 8. The dH/dt values over a sub-area of Jakobshavn Isbras, using average dH/dt values over the 
period 1999-2008. (a) Locations of ERS-2 and ICESat measurements and (b) interpolated ATM (consid- 
ered to be the “true” pattern). (c-f) Interpolation results by OK, KED, ST-OK, and ST-KED, in that order. 
Contours show velocity: dark and bright green represent 100 m yr -1 and 1000 m yr -1 velocity intervals, 
respectively. Only the areas further than 5 km from the grounding line and with velocities between 
100 m yr" 1 and 4000 m yr -1 are shown as these are used in most subsequent analyses (see text). 
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Table 3. Statistics of Data and Interpolations of dH/dt Values 
Over the Period 1999-2008 3 



Mean 

SD 

Min. 

10-p. 

90-p. 

Max. 

Input (760) 

-0.764 

1.115 

-6.989 

-2.157 

-0.008 

1.138 

ATM data (450) 

-2.413 

2.365 

-17.077 

-5.496 

-0.330 

0.565 

interpolated ATM 

-1.207 

1.877 

-16.688 

-3.041 

-0.112 

1.800 

OK 

-1.042 

1.356 

-6.823 

-3.265 

-0.091 

0.763 

KED 

-1.157 

1.776 

-19.863 

-3.513 

-0.081 

0.764 

ST-OK 

-0.731 

0.882 

-4.247 

-2.135 

-0.052 

0.721 

ST-KED 

-0.979 

1.793 

-26.124 

-2.638 

-0.052 

2.418 


a For space-time kriging, the value for this period was obtained by 
averaging dH/dt over individual years. Shown are the mean, standard 
deviations, minimum and maximum, as well as the 10% and 90% 
percentiles. Note that for input and ATM data, the numbers (between 
brackets) are very much smaller than the total number of pixels in the area 
of interest (14,348) and thus heavily influence the statistics. 

box-plots in Figure 9a. For spatiotemporal kriging, this 
exercise could only be carried out for individual years. For 
cross-validation only the input data (i.e., no ATM) were 
used. These were located almost exclusively outside the 
high-velocity area. The influence of the external drift should, 
therefore, be non-existent or very small. Figure 9a shows 
that this is indeed the case, except for some years prior to 
2002 where only (very few) ERS-2 measurements are 
available and one observation close to the high-velocity area 
has a noticeable effect. When comparing spatiotemporal 
kriging versus spatial-only kriging the differences are again 


small, although in most cases the errors are slightly smaller 
for ST-OK/ST-KED, because more (past and future) sam- 
ples are taken into account. Exceptions are 1998 and 2003. 
2003 is also the only year with a bias (i.e., the average of the 
errors is not zero). Because neither ERS-2 nor very many 
ICESat measurements are available for 2003 this can also be 
attributed to the low sampling density. 

[36] In Figure 9b, box plots are shown for a comparison of 
interpolated values and ATM measurements where those are 
available. Because the number of ATM measurements 
changes dramatically over the years, interpretation of these 
errors is difficult. To compare results for the period 1 999— 
2007, dH/dt from KED/ST-KED were averaged over these 
years. When many input measurements are available (typi- 
cally the middle of the ICESat era and in the case of the 
averaged data), using the external drift reduces both the bias 
(mean error) and the maximum errors. If only few input data 
are available, i.e., only ERS-2, either all methods produce 
nearly the same errors, or maximum errors are slightly larger 
when using the external drift. There are two cases where 
errors are larger for RED (2004) or ST-KED (2002). This 
period marks the end of the ERS-2 coverage and the 
beginning of ICESat coverage. Because 3-year running 
windows were used to select data for the regression, in this 
period very few dH/dt estimates with sufficiently low stan- 
dard error are available. Hence, the external drift “over- 
shoots" the ATM measurements in some of those occasions. 
Furthermore, it should be noted that in both years the 




Figure 9. (a) Errors of cross-validation and (b) comparison with ATM, for all individual years as well as 
the entire period (All). The numbers in Figure 9b represent the number of available ATM measurements. 
Furthermore, the size of the boxes represent the quartile range, with the middle being the median and 
the +-sign the mean. Whiskers indicate 1% and 99% percentiles, and the outer stars the minimum 
and maximum values. 
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Figure 10. Spatial pattern resulting from the four kriging methods (from left to right: OK, KED, ST-OK, 
ST-KED), for seven selected years (after 2005 the patterns are similar to 2005). The color scale shows 
dH/dt in m yr~ 1 . 


available ATM measurements were clustered in a small area 
(Figures 2c and 2e). 

[ 37 ] To further assess the temporal evolution of dH/dt, 
after interpolation. Figure 10 shows the difference in time 
evolution between the methods. Six selected years from the 
period 1997-2007 are shown. It should be noted that from 
2005 onwards the spatial patterns are similar. As mentioned 
before, prior to 2003 only data from ERS-2 was available, 
making it difficult to obtain a realistic pattern given the 
sparsity of ERS-2 (Figure 1). For KED and ST-KED, the sign 
of (sometimes very few) data points detennines the direction 
in which the trend is calculated. For example in 2003, where 
the few satellite measurements that are available are pre- 
dominantly positive (probably caused by a positive anomaly 
in the surface mass balance), values in the high-velocity area 
are predicted to be higher than the remaining part of the 
basin. This is not necessarily unrealistic, given the correlation 
between dynamic thickening and velocity that was discussed 
before (Section 3.3 and Table 1). 

[38] Using spatiotemporal kriging, where measurements 
from adjacent years are employed, smooths the temporal 
variations considerably. Figure 10 shows that the transition 
between thickening and thinning is around 2000, which is 
about two years later than reported elsewhere [e.g., Thomas 


et al, 2009; Csatho et al., 2008; Joughin et al., 2008b]. 
However, as mentioned before, only ERS-2 data were 
available in this period, none of which is located around the 
high-velocity area. The dynamic thinning started near the 
grounding zone and then propagated upstream. It thus took 
some time to advance far enough upstream to be captured by 
ERS-2. The ST-OK and ST-KED patterns around this tran- 
sition do not look realistic, as there is a sharp boundary 
between positive and negative values, caused by measure- 
ments with opposing signs. A higher sampling density in 
that period (e.g., additional sensors, or repeat-track proces- 
sing of ERS-2 instead of crossovers) would alleviate this. 


5. Discussion 

[ 39 ] In our results, kriging with external drift produces 
generally higher thinning rates that are closer to the ATM 
data. In addition, spatiotemporal kriging smooths effects of 
interannual variability in surface mass balance and it also 
improves estimates in periods with data gaps because it can 
use data from adjacent epochs. We believe, therefore, that 
ST-KED produces the best results. There are, however, some 
areas and periods where performance is less than optimal. 
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Table 4. Volume Change Estimates for Jakobshavn Isbrae"* 


Period 

OK 

KED 

ST-OK 

ST-KED 

1999-2007 

-13.6 

-15.2 

-7.6 

-11.1 

2003-2007 

-16.2 

-18.7 

-17.5 

-20.5 

1996 

18.1 

17.3 

16.2 

16.8 

1997 

9.0 

9.1 

9.0 

9.5 

1998 

-5.9 

-5.0 

-5.5 

-4.9 

1999 

5.9 

7.1 

4.7 

5.4 

2000 

6.1 

6.3 

5.2 

5.1 

2001 

5.0 

2.7 

1.6 

-5.3 

2002 

5.8 

2.7 

7.2 

-2.4 

2003 

-13.7 

-12.2 

-20.1 

-25.2 

2004 

-14.1 

-20.0 

-11.1 

-13.3 

2005 

-12.1 

-13.3 

-13.5 

-15.5 

2006 

-17.8 

-20.9 

-18.9 

-21.0 

2007 

-23.1 

-27.1 

-23.8 

-27.4 


a Volume change estimates are shown in km 3 yr 1 for individual years 
1996-2007, as well as averaged over the periods 1999-2007 and 2003-2007. 


[ 40 ] Close to the grounding zone, where velocity rapidly 
increases from ~4 km yr -1 to about ~12 km yr -1 in a dis- 
tance of 10 km, thinning rates as interpolated by KED/ 
ST-KED increase with a similar gradient, resulting in physi- 
cally unrealistic values. In this area, the relation between 
velocity and dH/dt breaks down because of uncertainty in 
the exact location of the grounding line. In addition, the 
large change in velocity gradient suggests that a different 
physical process becomes important. Furthermore, in the 
transition period from thickening to thinning, samples with 
opposing signs cause the trend from the external drift to 
have different signs, resulting in unrealistic patterns. This 
transition period, which in the high-velocity area happened 
around 1998, is delayed for some years in our interpolated 
time series because it took about two years for the thinning 
to propagate far enough upstream to be captured by radar 
altimetry. 

[ 41 ] Because KED and ST-KED yield higher thinning 
rates over the fast-flowing part of Jakobshavn Is bras, the 
resulting mass change might be different as well. The cal- 
culation of actual mass changes requires modeling density 
profiles and fim compaction, which is outside of the scope 
of this study. It is, however, straightforward to calculate 
volume changes by integrating dH/dt over the basin. Table 4 
provides volume changes for individual years between 1996 
and 2007, as well as longer-term averages. In these esti- 
mates, the area with velocity >4000 m yr - 1 and within 5 km 
of the 2004 grounding line was not taken into account, as 
discussed previously. These areas were assigned a dH/dt of 
—20 m yr -1 , which is a reasonable value given the minimum 
value in ATM of —24 m yr“ 1 . However, because this area is 
small, its impact on the integrated volume change is also 
small. As an extreme example, not assigning a value (i.e., 
0.0 m yK 1 ) yields volume losses that are about 2 km 3 yr^ 1 
lower. In general, volume change estimates for Jakobshavn 
Isbrae are about 10-20% higher for ST-KED than for ST-OK, 
with an average volume loss over 2003-2007 of 20.5 km 3 
yr -1 for ST-KED versus 17.5 km 3 yr“ 1 for ST-OK (Table 4). 
Thus, the interpolation technique used can have a significant 
impact on the estimated volume changes at this scale. How- 
ever, Jakobshavn Isbras is one of the largest glaciers in 
Greenland with high thinning rates over a large area. When 
the entire ice sheet is considered, the differences can be 


expected to be smaller because not all outlet glaciers are 
thinning. To assess the difference in volume or mass changes 
between interpolation methods for the entire ice sheet will be 
part of future work. 

[ 42 ] Furthermore, it would be useful to assess whether 
ST-KED is suitable for use on other outlet glaciers, or 
whether Jakobshavn Isbrae is particularly suited to this 
approach. An important requirement is that dynamically 
induced dH/dt values are dominant with respect to the SMB 
induced dH/dt. If that is not the case, it could be possible to 
use SMB from a regional climate model [e.g., Ettema et al, 
2009] to correct for the influence of SMB. If ice dynamics 
dominate, the assumption that significant dynamic thinning 
occurs mainly in the high-velocity area would allow the use 
of KED, as the exact relation between dH/dt and velocity 
(which is different for each glacier) is determined implicitly 
and is locally variable. In areas where spatial gradients in 
velocity are very low, the spatial trend in dH/dt will also be 
very small so that interpolated values will be very similar to 
those obtained by OK. Another requirement is that ice 
velocities are available at every location on the ice sheet. 
However, nearly all outlet glaciers are covered by InSAR 
measurements [e.g., Moon et al., 2012] and the gaps that 
occur are mainly located in the interior where gradients are 
small. 

6. Summary and Conclusions 

[ 43 ] Interpolation of dH/dt values derived from satellite 
altimeters is hampered by sparse sampling over, often, nar- 
row outlet glaciers, where the largest changes occur due to 
dynamic thinning. To mitigate this problem, we have used 
velocity data from InSAR and speckle tracking as secondary 
data to inform the interpolation. The approach was tested 
over Jakobshavn Is bra, where extensive airborne laser 
altimetry data ( ATM) are available for validation. 

[ 44 ] First, the relation between ice velocity and dH/dt was 
examined, revealing a near-linear relation with highly neg- 
ative correlation coefficients (<— 0.85). Interestingly, when 
the glacier was thickening in the mid-nineties a correlation 
was found as well, suggesting the thickening had a ice- 
dynamical origin (i.e., a slowdown at some earlier time). 
Using ATM data as reference, interpolated spatial patterns of 
average dH/dt are more realistic when velocity is used as a 
background field. Maximum thinning rates, that are about 
7 m yr” 1 for ordinary kriging, become similar to ATM data, 
at about 20 m yr -1 , when velocity is considered in the 
interpolation. In addition, spatiotemporal kriging smooths 
effects of interannual variability in surface mass balance and 
it also improves estimates in periods with data gaps because 
it can use data from adjacent periods. Errors with respect to 
ATM (mainly located in the high-velocity area) are generally 
reduced significantly, except for some cases with very few 
available observations. For these reasons, interpolation with 
ST-KED was considered to yield the best results. However, 
very close to the grounding line, spatial gradients of dH/dt 
become physically unrealistic due to a change in flow 
regime in this region. In addition, the transition from thick- 
ening to thinning (around 1998 in the high-velocity area) 
was delayed in our interpolated time series because it took 
time for the thinning to propagate into the areas with higher 
elevation where the ERS-2 radar altimeter was sampling. 
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[ 45 ] Volume changes that were calculated from the four 
resulting interpolation methods showed that volume loss 
derived from ST-KED was about 2-4 km 3 yr -1 , or 10-20%, 
higher compared to OK. When applied to the entire ice sheet 
ST-KED may therefore result in increased mass loss esti- 
mates. However, as Jakobshavn Is but' is one of the largest 
and most dynamically active glaciers in Greenland, the sig- 
nificance of this percentage increase will likely be smaller 
and remains for further research. 

[46] Previously, in order to estimate the Greenland mass 
balance, radar altimetry was merged with ATM data over 
outlet glaciers [Zwally et al., 2005]. The interpolation method 
proposed here yields similar results with radar altimetry only. 
Although for newer and future sensors such as CryoSat-Il 
and lCESat-2 coverage over outlet glaciers will likely be 
less of a problem, our interpolation will improve ice sheet 
mass balance reconstructions from existing and past satellite 
altimeter data sets. 
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